Ab-initio Molecular Dynamics 
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Computer simulations and molecular dynamics in particular, is a very powerful method to provide 
detailed and essentially exact informations of classical many-body problems. With the advent 
of ab-initio molecular dynamics, where the forces are computed on-the-fly by accurate electronic 
structure calculations, the scope of either method has been greatly extended. This new approach, 
which unifies Newton's and Schrodinger's equations, allows for complex simulations without relying 
on any adjustable parameter. This review is intended to outline the basic principles as well as 
a survey of the field. Beginning with the derivation of Born-Oppenheimer molecular dynamics, 
the Car-Parrinello method as well as novel hybrid scheme that unifies best of either approach are 
discussed. The predictive power is demonstrated by a series of applications ranging from insulators 
to semiconductors and even metals in condensed phases. 



INTRODUCTION 

The geometric increase in performance of computers 
over last few decades, together with advances in applied 
physics and mathematics has led to the birth of a new 
way of doing science that is in the intersection of the- 
ory as well as experiment. As a consequence, they are 
referred to as computational sciences and allow for com- 
puter experiments under perfectly controllable and re- 
producible conditions. In this way computer simulations 
have been very successful in explaining a large variety of 
physical phenomena and guiding experimental work. In 
addition it is even possible to predict new phenomena by 
conducting experiments in silico that otherwise would be 
too difficult, expensive, or simply impossible to perform. 
However, the by far most rewarding outcome of computer 
simulations is the invaluable insight they provide into the 
behavior and the dynamics of a system. The two most 
common algorithms for such studies are the Monte Carlo 
(MC) [H and Molecular Dynamics (MD) 0-0] algorithm. 
The latter is simply the numerical solution of Newton's 
equation of motion, which allows both equilibrium ther- 
modynamic and dynamical properties of a system at fi- 
nite temperature to be computed. Since it also provides 
a 'window' onto the atomic real-time evolution of the 
atoms, another role of MD is that of a computational 
microscope. 

One of the most challenging, but very important as- 
pects of MD simulations is the calculation of the inter- 
atomic forces. In classical simulations, they are com- 
puted from empirical potential functions, which have 
been parametrized to reproduce experimental or accu- 
rate ab-initio data on small model systems. Even though 
great strides in elaborating these empirical potentials 
have been made, often the transferablility to systems 
or regions of the phase diagram different from the ones 
to which they have been fitted is restricted. Further- 
more, they are not able to simulate with sufficient pre- 
dictive power chemical bonding processes that take place 



in many relevant systems. Eventually, some of the most 
important and interesting phenomena of modern physics 
and chemistry are intrinsically nonclassical. 

Therefore, a first-principles based approach, such as 
ab-initio molecular dynamics (AIMD) |5(, where the 
forces are calculated on-the-fly from accurate electronic 
structure calculations, is very attractive since many of 
these limitations can in principle be removed. However, 
the increased accuracy and predictive power of AIMD 
simulations comes at significant computational cost. For 
this reason, density functional theory (DFT) is to date 
the by far most commonly employed electronic structure 
theory, but it is important to note that AIMD is a general 
concept that in principle can be used in conjunction with 
any electronic structure method. Nevertheless, also the 
ab-initio approach is not without problems - the relevant 
energy scale is tiny, well below fcgT, and in particular 
the attainable length and time scales are still one of its 
major limitations. 



MOLECULAR DYNAMICS 

The mathematical task of MD is to evaluate the expec- 
tation value (O) of an arbitrary operator 0(R,P) with 
respect to the configurational Boltzmann distribution 



(O) 



JdRdPOjR, p) e -mn,P) 
JdRdPe-P E ( R : p ) 



(1) 



where /3 = is the inverse temperature. The 

total energy function 



E(RP) 



N io „ 



Pf 



^ 2Mi 

i=i 1 



(2) 



where the first term denotes the nuclear kinetic energy, 
Q(Ri) the potential energy function, Ni on the number of 
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ions and Mj the corresponding masses, depends itself on 
nuclear positions R and momenta P. 

One way to to evaluate Eq. ([I} , at least in principle, is 
to directly solve such a high dimensional integral, whose 
integrand is very sharply peaked in many dimensions, by 
an uniform sampling using the MC algorithm. However, 
such an MC algorithm is very inefficient, if it would not 
be for importance sampling 1], which satisfies the suffi- 
cient detailed balance condition by rejections. 

On the other hand, assuming the ergodicity hypothe- 
sis, the thermal average (0) can not only be determined 
as the ensemble average of a MC simulation, but using 
MD equally as a temporal average 

(0)= Hm i f dtO(R(t),P(t)). (3) 

T->oo / J 

However, by propagating the classical many-body sys- 
tem in time according to Newton's equation of motion, 
the ions are treated only classical, an usually negligible 
approximation except for very light atoms or low temper- 
ature, where nuclear quantum effects may be important 
and a quantum formalism such as imaginary-time path 
integrals 0, Sj is required. 

Similar to MC, within MD some kind of importance 
sampling is naturally performed by preferentially visit- 
ing phase space of low potential energy. Furthermore, as 
denoted by the additional time dependence in Eq. (|3|), 
MD allows for additional insights from the ionic real- 
time evolution, at least in an statistical average sense. It 
is neither the intention, nor possible, to obtain "exact" 
trajectories by MD due to the infamous Lyapunov insta- 
bility, which states that slightly perturbed trajectories 
are intrinsically exponentially diverging with time. 

The equipartition theorem 

(-MtR 2 ^ = h B T, (4) 

where Mj is the atomic mass, ks the Boltzmann con- 
stant and T the instantaneous temperature, offers an ele- 
gant way to bridge the gap between molecular mechanics 
and thermodynamics. This opens the door to extract a 
vast variety of relevant static and dynamic, as well as 
transport properties from a MD simulation. 

Nevertheless, any computational resource is finite, 
which limits the time and length scales accessible by 
computer simulations. One way to partially bridge the 
gap between the microscopic size of the simulated sys- 
tem and the macroscopic reality is to introduces periodic 
boundary conditions. In this way surface effects are elim- 
inated by effectively simulating an infinite system, albeit 
with a finite periodicity that is identical with the length 
L of the simulation cell. As a consequence, only phe- 
nomena whose characteristic correlation length is much 



smaller than L can be simulated. By similar means only 
processes whose typical relaxation time is significantly 
smaller than the simulation time T can be studied. Even 
though great strides have been made to extend the time 
and length scales, it is apparent that techniques such as 
those reviewed here are clearly needed. 

AN AB-INITIO POTENTIAL 

In AIMD the forces Fj = — Vi^ <&(!?/) are determined 
on-the-fly using first principles electronic structure meth- 
ods. That means that AIMD is not relying on any ad- 
justable parameter, but only on Rj, which constitutes 
its predictive power. However, finding the antisymmet- 
ric ground state eigenfunctions ipo of the corresponding 
many-body Hamiltonian at each MD step comes at a 
significant computational cost, which has to be carefully 
balanced against the size and sampling requirements of 
MD. 

The Many-Body Schrodinger Equation 

Applying the so called Born-Oppenhcimcr approxima- 
tion [9( , which we have implicitly assumed in the preced- 
ing subsection, $(_Rj) can written as 

= (ih\Ue({ri}; RiMo) + Ejj(Rj), (5) 

where H e ({ri};Rj) is the electronic many-body 
Hamiltonian, that depend on the electronic coordinates 
{r.;} and parametrically on R]. Essentially, the Born- 
Oppenheimer approximation allows for a product ansatz 
of the total wavefunction consisting of the nuclear and 
electronic wavef unctions. Due to the large separation 
of the nuclear and electronic masses, the electrons can 
be expected to be in its instantaneous equilibrium with 
the much heavier nuclei, so that the electronic subsys- 
tem can be treated independently at constant R]. Nev- 
ertheless, we are left with a formidable task to solve the 
electronic, nonrelativistic, time- independent, many-body 
Schrodinger equation 

% e {{n}; Ri)M{n}) = eo(Ri)M{n})> ( 6 ) 

which is a high-dimensional, non-linear eigenvalue 
problem, with eigenfunctions T/>o({ r i}) and eigenvalues 
Eq(Ri), respectively. 

To visualize the complexity of Eq. (|6]) let us consider 
the following Gedankenmodell to represent the solution 
ipo({ r i}) ° n a real-space grid, where each coordinate is 
discretized by as less as 100 mesh points. Ignoring spin 
and taking ipo({ r i}) to be real, for N e electrons 10 6JVc 
grid points are required, so that the solution of a single 
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Si atom would require more grid points than the number For the sake of simplicity in the following I will 
of electrons in the whole universe, not to mention solving throughout assume atomic units and confine myself to 
such a large non-linear eigenvalue problem. the physical relevant coulomb system, for which 



Density Functional Theory 

Fortunately, this curse of dimensionality can be inge- 
niously bypassed by the use of DFT, which is based on 
two celebrated papers of Hohenberg, Kohn and Sham 
6, 10]. The former, the so called Hohenberg-Kohn (HK) 
theorem, proves the existence of an one-to-one mapping 
between the ground state density po( r ) and an external 
potential v(r). In this vain, po(r), which depends on just 
3 electronic degrees of freedom, is designated as the prin- 
cipal quantity rather than the 3-/V e -dimensional many- 
body wavefunction. As a consequence, the nondegener- 
ate ground state wavefunction ipo({ri}) — ip[po(r)] and 
likewise H e [po( r )] are both unique Junctionals of po{r), 
just as the ground state energy Eo = E BFT [po(r)] = 
(tp{po(r)}\'H e [po(r)}\ip{po(r)]) . The latter obeys the vari- 
ational property 
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Po] = (^o|H e |^ ) < W\U e W) = E \p'\ i (7) 



for which equality holds if and only if po = p'. As a 
consequence Eq. ((6]) can not only be solved by iteratively 
diagonalizing H e [p] within a self-consistent field (SCF) 
procedure, but equally by minimizing the quantum ex- 
pectation value of H e [p\- 



E DFT [p ] = min(V>|H e |^) = mm ty\p]\He[pM\p]) (8a) 

= min J B DFT [ j o]. (8b) 
p 

In principle the minimization have to be performed 
under the constraint that p(r) is ./V-representable, i.e. 
that it is arising from an antisymmetric N-body wave- 
function ip({ri}). Luckily, this had been solved, and 
it can be demonstrated that any single-particle den- 
sity can be written in terms of an antisymmetric many- 
body wavefunction 
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Il2j . On the contrary, for the 
w-representability problem, which states that p(r) is the 
ground state density of a local potential v(r), no such 
general solution is known. The HK theorem just guaran- 
tees that there cannot be more than one potential for each 
density, but does not exclude the possibility that there 
is no potential realizing that density. It is only known 
for discretized systems that every density is interacting 
ensemble w-representable. Interestingly, the constructive 
proof of Levy and Lieb [HI, EH shows that for an inter- 
acting system i>-representability is not required for the 
proof of the HK theorem. 



i—1 i<Cj 

= f + u + v, 



(9b) 



where T is the kinetic energy operator of the elec- 
trons, while U is the electron-electron interaction and 
V = vfri) the electron-ion operator. The former two 
operators are universal and independent of the system, 
while the latter is system dependent, or nonuniversal. 
DFT explicitly recognizes that it is indeed the potential 
v(r), which distinguishes nonrelativistic Coulomb sys- 
tems and offers a prescription how to deal with T and 
U once and for all (Isj . Hence, even at this stage based 
using nothing but the HK theorem DFT is already of 
some practical use without having to solve the many- 
body Schrodinger equation and without having to make 
a single-particle approximation. In principle it should be 
even possible to calculate all observables, since the HK 
theorem guarantees that they are all functionals of p(r). 
Presuming the availability of physical sound and accurate 
approximations one can write 
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[p(r)} = T[p(r)} + U[p(r)} + V[p(r)}. (10) 



In the so called Thomas-Fermi (TF) approximation 
fl6l 17 ] the full electron-electron interaction energy is 
approximated by the Hartree energy 



U[ P {r)]^U H [p{r)]= l - j dr J dr 



, P(r)p(0 (u) 
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that is the electrostatic interaction energy of p(r). In 
addition the kinetic energy is approximated as 



T\p(r)] 



drt nom (p(r))=T LDA [p(r)} 7 (12) 



where t hom (p(r)) is the kinetic-energy density of a ho- 
mogeneous interacting system, which is also known as 
the local-density approximation (LDA). Due to the fact 
that the explicit form of t hom (p(r)) is only known for a 
non- interacting system, t hom (p(r)) is further estimated 
by the single-particle approximation t^ om (p(r)), i.e. 



T LDA [p(r)]« J 'drt h s om (p(r))=T^ A [p(r)}. (13) 
In the end the TF energy functional 
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£ TF [ P (r)] = r s LDA [p(r)] + UhW)\ + VW)\ (14) 

implies not only the single-particle approximation to 
the full electron-electron interaction, but also the single- 
particle mean- field approximation T s LDA [p(r)] to the ex- 
act kinetic energy of the inhomogeneous interacting sys- 
tem. As a consequence, all many-body correlation effects 
are neglected. 

However, the HK theorem predicates that they arc 
again a functional of p(r). The addition of an approxima- 
tion to the exact exchange and correlation (XC) energy 
results in a formally exact theory, which is referred to as 
orbital- free DFT [18j]. It is therefore important to rec- 
ognize, that the HK theorem is nothing but the formal 
exactification of the TF approximation. 

Similarly, the Kohn-Sham (KS) scheme can be con- 
sidered as the exactification of the self-consistent Hartree 
equations (HE) [19], which differs only in the kinetic en- 
ergy from the TF approximation. In fact, for the ficti- 
tious non-interacting system the kinetic energy is known 
exactly, even though only in terms of an explicit single- 
particle orbital functional, i.e. as an implicit density 
functional 



Therewith, the KS energy functional is simply given 

by 

E KS [p(r)} = E KS [WMr)]}] = T s [{^\p(r)]}] 

+ U H [p{r)\ + V[p(r)} + E xc [p(r)} (18a) 

N 

2 J \r-r'\ 
+ J drv^ t (r)p(r) + E xc [p(r)}, (18b) 

where E xc [p(r)} = (T[p(r)] - T s [{i/Ji[p{r)}}}) + 
(U[p(r)] — U~H[p{r)]) is the already mentioned and ap- 
parently unknown XC energy functional, whereas v X c ( r ) 
is the corresponding XC potential. This definition also 
shows that a significant of part Exc[p(r)] is due to cor- 
relation effects of the kinetic energy, that is known ex- 
plicitly only in terms of the reduced 2-particle density 
matrix [2pj | . 

Since up to the exactifying term E xc [p(r)} Eq. (|18a[) 
is identical to the HE, it is not surprising that the corre- 
sponding Euler-Lagrangian equation 



T s [p(r)} =-oE / <*r#(r)VV<(r) (15a) 

= T S [{^ [/>(*■)]}]■ (15b) 

Here the fictitious single-particle orbitals, or simply KS 
orbitals, are denoted as ipi(r). As we will see immediately 
they are eigenfunctions of a fictitious system, known as 
the KS system. It therefore should be noted that they dif- 
fer from the single-particle orbitals used in wavefunction 
based methods and have no strict physical meaning, with 
two notable exceptions: at the presence of the exact XC 
functional for the special case of an isolated system with 
v(oo) = (i) the highest occupied eigenvalue Sn can be 
shown to be the negative of the exact, many-body, first 
ionization potential including relaxation effects, and (ii) 
that the lowest unoccupied eigenvalue e^+i is the nega- 
tive of the electron affinity. Beside these two exceptions, 
only the density has a real physical meaning and can be 
written in terms of ipi( r ) as 



N 



(16) 



i=l 



where N occ is the number of occupied orbitals and fi 
the occupation number of state i, so that 



N 



(17) 



--V 2 +T, s KS (r))<Mr) 



(19) 



also results in a similar fictitious single-particle equa- 
tion. Since v^ s (r) — vn(r)+v xc (r)+v(r) is the effective 
potential of an artificial system, such that the ground 
state density and therewith the energy equals those of 
the true interacting many-body system. This particular 
system is therefore called KS system, its effective poten- 
tial KS potential and the resulting set of self-consistent 
equations are referred to as KS equations: 



-V 2 + vg.(r) + vxc(r) + v(r) ipi(r) = e l ip i {r) (20a) 



N„ 



Y,Mi(rW*(r) = p(r) (20b) 



SE xc [p(r)} 
5p{r) 



v xc (r) (20c) 



At self-consistency it is possible to express E^ s [p(r)] in 
terms of the single-particle KS eigenvalues Due to the 
fact that they are not the eigenvalues of the interacting 
many-body system, but of fictitious non-interacting KS 
system, E^ s [p(r)] is merely the sum of Si, but 



E$ S [p(r)] = J2fei-2 J drdr 



, P(r)p{r') 



\r — r 



J drv xc (r)p(r) + E xc [p(r)}. (21) 
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That is to say that in order to make genuine calcula- 
tions the KS scheme systematically maps the full inter- 
acting many-body problem, with U, onto an equivalent 
fictitious single-body problem, with an effective potential 
operator Vks — U s + Vh + Vxc, but without U [lBj : 



TF Exc[pir) \ HK 



r 8 [«dpO)]}]| 
HE 



Exc[p(r)] 



Ts[Wi[p(r)]}] 

-> KS 



The Exchange and Correlation Functional 

In the previous subsection DFT has been outlined as 
an exact theory, presuming that the exact XC functional 
is known. Unfortunately, except for the uniform electron 
gas , this is not the case and one has to resort to more 
or less accurate approximations. For the sake of brevity, 
only the main physical principles, as well as the most 
notable exact properties, rather than the various rungs 
of " Jacob's ladder to heaven" [22j will be discussed here. 

On this account the following break-up is particularly 
convenient: 



Exc [p(r)\ = i J dr J dr 



, p(r)pxc(r,r') 



\r — r' 
E x [p(r)] + E c [p(r)}, 



(22) 



where Ex[p(v)] is the exchange energy due to Pauli 
repulsion, Eq [p(r)] the electron correlation energy and 
Pxc{ r , r ') = Px{r 7 r') + pc{r 7 r') the XC hole. The for- 
mer is therefore the energy lowering due to the antisym- 
metry requirement on the wavefunction of a fermionic 
system and can be exactly calculated in terms of an ex- 
plicit orbital, or implicit density functional 



Ex[p(r)] = -§ > / «''• 

hi 

= E x [mp(r)}} 



dr 



,^(r)^(r)Vi(r)^(r)' 



(23) 



Eq. (|23|) is the so called Hartree-Fock exchange energy, 
but with KS orbitals. However, the nonlocal form of the 
exact exchange energy comes at a considerable compu- 
tational burden to solve four-center integrals, which is 
about two order of magnitude more expensive than is 
the case for local or semi-local approximations to the ex- 
act XC functional. The correlation energy accounts for 
the additional energy lowering, since electrons with op- 
posite spins also avoid each other. However, contrary to 
the exchange part, no exact expression for Ec[p{r)] is 
known, neither in terms of orbitals nor densities. 



Obviously, DFT would be of little use if one had to 
know Exc[p{ r )] exactly, but luckily it is usually ener- 
getically substantially smaller than each the remaining 
terms, which are known. One can thus hope that reason- 
able simple approximations to Exc[p( r )], will still allows 
for qualitatively correct estimates of E [p(r)], without 
relying on additional adjustable parameters. 

The following properties of the exact XC functional 
are well established and serves as valuable guides in 
the construction of approximations to Exc[p(f)]- In 
particular the exact exchange hole px(r,r') obeys the 
sum rule J dr' px(r,r') = —1, and since the correla- 
tion hole integrates to zero furthermore / dr' px{r, r') — 
J dr' pxc( r : r ') is valid. 

Another important property of the exact XC functional 
is that in the one-electron limit the following conditions 
hold: 



Ec[p (1) (r)} = (24a) 
E x [p W (r)]=-E H [pW(r)}, (24b) 

where p^ is a one-electron density. These conditions, 
which are satisfied within the Hartree-Fock approxima- 
tion, but not by standard local or semilocal XC function- 
als, ensure that there is no artificial self-interaction of an 
electron with itself. 

The Lieb-Oxford bound [H 



ExW)] > ExcW)\ > -1-68 J rfrp(r) 4 / 3 , (25) 

represents a lower bound on the XC energy and is typ- 
ically fulfilled by the most common XC functionals. 

One, if not the most, intriguing property of the exact 
XC functional, is the derivative discontinuity of the XC 
functional with respect to the particle number 24T 26j| 



Axe = 



5ExcW)] 
5p{r) 



6Exc[p{r)} 
5p{r) 



J XC 



(r) 



J xc 



(26) 



where 5 is an infinitesimal shift of N ei while Axe is a 
system-dependent shift of vxc( r )j which is independent 
on r. The noninteracting, single-particle kinetic energy 
functional obeys a similar discontinuity given by 



ST s [mp(r)}}} 
5p{r) 

£N+i{r) - s N (r), 



ST s [mp(r)}}} 
5p{r) 



N c -S 

(27) 



where En+i and are again the KS single-particle 
energies of the lowest unoccupied and highest occu- 
pied eigenstate. In the quantum chemistry community 
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they are denoted as lowest unoccupied molecular or- 
bital (LUMO) and highest occupied molecular orbital 
(HOMO), respectively. Therefore, Eq. (J2ZJ) is nothing 
but the KS single-particle gap A K s, or HOMO-LUMO 
gap, whereas Eq. (|26p is a truly many-body effect. There- 
fore, the true fundamental gap A = E(N+1) - 2E(N) + 
E(N-l) is the derivative discontinuity of the exact total 
energy functional for the ground state [2(| 



A 



5E[p(r)} 
= Aks + Axe- 



SE[p(r)] 
5p(r) 



(28) 



Since all other terms in E[p(r)] are smooth functionals 
of p(r), their derivatives are continous, so that A is only 
the sum of Eqs. (|26[) and (f27|) . However, in particular 
the derivative discontinuity of the XC functional has till 
now resisted all attempts of describing it within local or 
scmilocal approximations. In fact, standard XC function- 
als predict Axe = 0, which causes that the fundamental 
gap is typically underestimated. 



AB-INITIO MOLECULAR DYNAMICS 

In the following let us assume that the potential en- 
ergy function is calculated on-the-fly using DFT, so that 
*(Jlj) = E[{^ i } l m]=E^ s [{ij i [p(r)}}]+Eu(Ri). In 
any case, AIMD |5j, l27H31j comes in two fundamental 
flavors, which are outlined in this section. 



one obtains the associated equations of motion (EOM) 



MiRj = -V Rl 
dE 



dRr 



nanE[{fl>i};Ri] 



{(^|V 3 >=««}. 



1,3 



d(il>i 



dRi 



SE 



(31a) 



0< 



SE 



■H e {iii\ 



(31b) 



The first term on the righ-hand side (RHS) of Eq. (|3"Ia|) 
is the so called Hellmann-Feynmanforce. The second 
term, which is denoted as Pulay (32[, or wavefunction 
force -Fwf , is a constraint force due to the holonomic or- 
thonormality constraint, and is nonvanishing if and only 
if the basis functions <j>j explicitly depends on Rj. The 
final term stems from the fact that, independent of the 
particular basis set, there is always an implicit depen- 
dence on the atomic positions through the expansion co- 
efficient dj (r) within the common linear combination of 
atomic orbitals <A,: 



(r) 



(32) 



Born-Oppenheimer Molecular Dynamics 

In Born-Oppenheimer MD (BOMD) the potential en- 
ergy E [{ipi}; Ri] is minimized at every MD step with 
respect to {ipi(r)} under the holonomic orthonormality 
constraint (ipi(r) \ Tpj{r)) = Sij. This leads to the follow- 
ing Lagrangian 



1 N 

C BO ({ipiY, Ri, Ri) = o MrRj - min£[{^}; Rj] 

1 = 1 XYtJ 

+ ^A y ((^ | V^)-%), (29) 



where A is a Hcrmitian Lagrangian multiplier matrix. 
By solving the according Euler-Lagrangian equations 



d dC 



dC 



dt dRj dRi 
d dC dC 



(30a) 
(30b) 



The factor 2 stems from the assumption that the KS 
orbitals are real, an inessential simplification. Never- 
theless, the whole term vanishes whenever ipi(r) is an 
eigenfunction of the Hamiltonian within the subs pac e 
spanned by the not necessarily complete basis set 3^, 34 1 . 
Note, that this is a much weaker condition than the 
original Hellmann-Feynman theorem 35, 36|, which we 
hence haven't availed throughout the derivation, except 
as an eponym for the first RHS term of Eq. (I31a[) . How- 
ever, as the KS functional is nonlinear, eigenfunctions 
of its Hamiltonian H e are only obtained at exact self- 
consistency, which is why the last term of Eq. (|31al) is 
also referred to as non-self-consistent force -Fnsc- Un- 
fortunately, in any numerical calculation this can not be 
assumed, which results in immanent inconsistent forces 
and to the inequality of Eq. (I31b[) . Neglecting either Fwf 
or -Fnsc j i-e. applying the Hellmann-Feynman theorem 
to a non-eigenfunction leads merely to a perturbative es- 
timate of the generalized forces [37[ 



F — -Fhf + i*WF + -Fnsc, 



(33) 
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which, contrary to the energies, depends just linearly 
on the error in the electronic charge density. As a con- 
sequence it is much more exacting to calculate accurate 
forces than total energies. 

However, as a corollary of the BO approximation, the 
electronic, as well as the ionic subsystems are adiabat- 
ically strictly separated from each other, and therefore 
does not entail any restrictions on the maximum possi- 
ble integration time step, so that time steps up to the 
ionic resonance limit are feasible. This actually holds ir- 
respective of the band gap, so as to, at least in principle, 
even metals can be straightforwardly treated. 



Car-Parrinello Molecular Dynamics 

In Car-Parrinello MD (CPMD) ^ a coupled electron- 
ion dynamics is performed, in which the electronic de- 
grees of freedom are added to the Lagrangian as classical 
ones: 



M N 



1=1 



E[{^};Ri] 



(34) 



Once again, applying the Euler-Lagrangian equations 
Eq. ([30af and Eq. (|30bj) entails an EOM, where the elec- 
tronic degrees of freedom inhere an artificial inertia /i and 
are propagated within a fictitious Newtonian dynamics, 
such that the electrons follows the ions adiabatically: 



MjRi 



dE 
SE 



d 



H e (i/3i\ + ^2A ij \ip j 



{(V> 4 hM=««}. 



(35a) 



(35b) 



As a consequence of the BO approximation, the high- 
frequency oscillations of Eq. (|35bl) vanishes on ionic 
time scales, so that ipi ~ 0. Hence, similar to Ehrcn- 
fest dynamics [38|, the total derivative of the instanta- 
neous, rather than the fully minimized, expectation value 
(^ol-ffel^o) °f the Hamiltonian yields the forces that are 
consistent with the corresponding energies. This means, 
that owing to the absence of necessity to fully minimize 
the energy functional, but rather to simply evaluate it at 
the instantaneous time step, Fnsc is identical zero by its 



very definition. Given a sufficiently small fictitious mass, 
the constant of motion is strictly conserved and errors in 
the forces are negligible, in particular if the ionic masses 
are renormalized by a constant mass tensor [HEI. In 
this respect CPMD combines most of the vantages of BO 
and Ehrenfest MD in the sense that the KS functional 
are only evaluated. There is no need to repeatedly solv- 
ing it cither by diagonalization or, equivalently, iterative 
minimization. However, due to the finite accuracy of any 
integrator, the holonomic orthonormality constraint of 
the orbitals have to be explicitly enforced. 

In order to ensure an adiabatic energy-scale separa- 
tion of the nuclear and the electronic degrees of freedom 
and to prevent energy transfer between them, the ionic 
phonon frequency Wj has to be much smaller than the 
electronic analog uj e . As had been suggested in an emi- 
nent phenomenological study of Pastore and Buda [42| 



AE, 



9«P 



(36) 



As a consequence the maximum integration time step 
Atmax depends on the inertia like ^JJi. The same also 
holds for the deviation from the BO surface 43 1 



|</v(r,t)-Vo(r,i)| <Cy/ji. 



(37) 



Therefore the fictitious mass, although physically com- 
pletely meaningless, acts as a continuous slider which al- 
lows to adjust any desired degree of accuracy, in terms 
of deviation from the BO surface, reciprocal to the com- 
putational efficiency in a well controlled manner. 

But if a metallic system is treated, due to the fact 
that CP states are strictly not KS eigenstates, Eq. ([36]) 
is identical zero and either a thermostat for the electronic 
degrees of freedom 3il 44, 4|| to counterbalance the ex- 
change of energy, or an extended functional with frac- 
tional occupation numbers [46j-|49( is necessary. 

In the end drawing a proper conclusion, if either 
BOMD or CPMD is to favour, turns out to be very sub- 
tle (Ho| and depends largely on the definition of accuracy 
and on the particular application. 



AN EFFICIENT AND ACCURATE 
CAR-PARRINELLO-LIKE APPROACH TO 
BORN-OPPENHEIMER MD 

Even though DFT-based AIMD has been very success- 
ful in describing a large variety of physical phenomena, 
its computational cost has limited the attainable length 
and time scales in spite of substantial progress. For a 
while it was believed that linear scaling methods [Sir- 
53j could have offered a solution. Unfortunately, the 
crossover point at which linear scaling methods become 
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advantageous has remained fairly large, especially if high 
accuracy is needed 



54, 



55| . Therefore, it would be very 
desirable to accelerate ab-initio simulations with up to 
thousands of atoms, such that simulations as long as a 
few nanoseconds can be routinely performed, thus mak- 
ing completely new phenomena accessible to AIMD sim- 
ulations. BOMD, in which the DFT functional is fully 
minimized at each MD time step, does not seem to of- 
fer much room for further improvement. For this reason 
recently another direction has been followed to improve 
the efficiency at current system sizes. In the spirit of 
CPMD [27|, some form of dynamics for the electronic 
degrees of freedom is implemented, which automatically 
keeps the system close to the instantaneous BO surface, 
but at variance to the original proposal in a localized or- 
bital representation 56|-59j . The acceleration stems on 



the one hand from this more compact description of the 
electronic wavefunctions, but is mainly due to the ability 
to reduce or even fully bypass the aforementioned SCF 
cycle. Nevertheless, just like in CPMD, these methods 
suffer from rather short integration time steps. However, 
rather recently a novel Car-Parrinello-like approach to 
BOMD has been proposed, which overcomes this limita- 
tion and combines the accuracy and long time steps of 
BOMD with the efficiency of CPMD [f^. 

From now on the general case will be considered, where 
the DFT KS orbitals are expanded in a non-orthogonal 
basis set. Let M be the dimension of the Hilbert space, 
i.e. the number of basis functions, and S the M x M 
overlap matrix. As usual the expansion coefficients of the 
N lowest occupied orbitals are arranged in a rectangular 
M xN matrix C. The density matrix can then be written 
as P = CC T and must obey the idempotency condition 
P = PSP that is due to the fermionic nature of electrons, 
which compels the wavefunction to be antisymmetric in 
order to meet the Pauli exclusion principle. The potential 
energy surface on which the ions move is defined by the 
minimum of an appropriately chosen energy functional 
-Edft [C, Rj], which is expressed as a functional of C and 
a function of the ionic coordinates Rj . In this notation 
the BO EOM reads as follows: 



posed. In Car-Parrinello-like approaches this is circum- 
vented by the design of a coupled electron-ion dynamics, 
which maintains the system very close to the BO surface, 
but at the cost of small integration time steps. 



Density Matrix Propagation 

Based on ideas of the original CP approach it is pos- 
sible to design an improved dynamics for the coupled 
system of electrons and ions [60(. However, contrary to 
the original scheme, this novel method is not expressed 
as an explicit EOM for the C's, but rather as an integra- 
tion scheme for the electronic degrees of freedom. The 
knowledge of the previous K values of C (t n -z), where 
I € [lf-K - ], determines the value of C(t n ), such that at 
any instant of time the C's are as close as possible to 
the instantaneous ground state. As for the short-term 
integration of the electronic degrees of freedom, accu- 
racy is crucial, a highly accurate and efficient algorithm 
is required. Therefore, here the always stable predictor- 
corrector (ASPC) method [sj EH of Kolafa has been 
selected. This scheme was originally devised to deal with 
classical polarization, so that care must be taken that 
during the evolution the idempotency condition is always 
satisfied. The modified predictor 



K I IK \ 

C p (t n ) = ^(-l) m+1 m)|^-C(t„_ m )C T (i n _ m ) 

{ K-l ) 



m— 1 

xS(t 

n—m )C(tn_i) 



p(t„. 



(39) 



uses the extrapolated contra-covariant density matrix 
PS as an approximate projector on to the occupied sub- 
space C(i„_i). In this way, the fact that the physi- 
cally relevant contra-covariant density matrix PS evolves 
much more smoothly and is therefore substantially eas- 
ier to predict than C is ideally utilized. The modified 
predictor is followed by a corrector step to minimize the 
error and to further reduce the deviation from the instan- 
taneous ground state. The corrector 



MjRi 



- V/ min -Edft [C, R/1 
c 



(38) 



where the search for the minimum is restricted to the 
C's that satisfy the orthonormality condition C T SC = I, 
which is equivalent to imposing the idempotency con- 
dition on P. As before, the forces of Eq. ([38]) can 
be divided into three contributions, (i) the Hellmann- 
Feynman forces 35, 36|, (ii) the Pulay forces 32 1, which 
are present whenever the basis set depends on the ionic 
positions, and (iii) a residual term [37[, which is non- 
zero except when full self-consistency is reached. The 
last term leads inevitable to a poor energy conservation 
in BOMD unless a very tight convergence criterion is im- 



C (*„) = wMIN [C p (t n )} + (1 - u)C p (t n ) 
K 

with 



2K- 1 



(40) 



consists only of a single preconditioned minimization 
step MIN [C p (<;„)] of a properly selected minimization 
procedure. Apparently, the predictor can also be repeat- 
edly applied, in which case the ground state is even more 
closely approached, but at the cost of additional elec- 
tronic gradient calculations. However, as will be shown 
immediately in general this is not necessary. The numer- 
ical coefficients of Eq. (|3T))) were selected in order to en- 
sure time-reversibility up to 0(h K+2 ), while u> was chosen 
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to guarantee a stable relaxation towards the minimum. 
Due to the fact that the energy is invariant under unitary 
transformations within the subspace of occupied orbitals 
C, it must be ensured that this time-dependent gauge 
transformation is not strongly changed by MIN [C p (t n )], 
as in this case continuity between the C's may be lost. 



Electronic Forces by Orbital Transformations 

Moreover, the minimization scheme must be very effi- 
cient in bringing the system as close as possible to the in- 
stantaneous ground state and at the same time preserves 
the idempotency condition of the density matrix. For 
these reasons, the orbital transformation (OT) method 



of VandeVondele and Hutter [63j has been chosen. In- 
spired by the form of the exponential transformation [64[ 
an auxiliary variable X is introduced, to parameterize 
the occupied orbitals 



C (X) = C p (t n ) cos (U) + XLT 1 sin (U) , 



(41) 



1 /2 

where U = (X T SX) and the variable X has to obey 
the linear constraint X T SC P (t n ) = 0. Under this condi- 
tion C (X) leads to an idempotent density matrix for any 
choice of X, provided that the reference orbitals C p (t n ) 
are orthonormal. Thus, any finite step along the precon- 
ditioned gradient direction will exactly fulfill the idempo- 
tency constraint by construction. Due to the linear con- 
straint the minimization with respect to X is performed 
in an auxiliary tangent space. Since this space is linear, 
no curved geodesies must be followed, as is the case for 
variables such as C that are nonlinearly constrained. In 
this way, large minimization steps can be taken, espe- 
cially if a good preconditioner is used 65|. In fact, us- 
ing an efficient, idempotency conserving direct minimizer 
such as OT is decisive for the success of this approach. 
Since the ASPC integrator only approximately preserves 
the idempotency constraint, it sporadical has to be ex- 
plicitly enforced, either by a Cholesky decomposition or 
by single purification iterations [661 ] . 



Total Energies and Forces 

Having obtained the new wavefunction it is now possi- 
ble to evaluate the energy and the nuclear forces, which 
are derived from the following approximate energy func- 
tional: 



E PC [p p ] = Tr [C T H[p p ] C]- l -JdvJ dr 
- J drV xc [p P ]p p + E xc [p p ] + E n , 



, p p (v)p p (r>) 



r — r 



(42) 



where p p (r) is the density associated with C p (t n ). 
Epc[p] can be thought of as an approximation to the 



Harris- Foulkes functional 67, \6&i and maintains the 



predictor-corrector flavor of this method. The validity of 
Epclfi] depends only on the efficiency of the minimizer 
and on the quality of the propagation scheme. The ionic 
forces are calculated by evaluating the analytic gradient 
of Epc[p] with respect to the nuclear coordinates. How- 
ever, as Ap = p — pP 7^ 0, besides the usual Hellmann- 
Feynman and Pulay forces an extra term appears: 



dr 



dVxc[p p ] 
dpP 



Ap + V H [Ap] 



(Vjp*H, (43) 



where p is the corrected density evaluated using C(t n ) 
and ff is the predicted density calculated from C p (t n ). 
Using variational density functional perturbation theory 

silza, E q . @3) can be efficiently computed very similar 

to employing the coupled-perturbed KS scheme. How- 
ever, due to the fact that usually only a single precondi- 
tioned minimization step is performed, C(t„) is just an 
approximate eigenfunction of H[p p ] within the subspace 
spanned by the finite basis set used. This leads to an 
insignificant error in the forces, provided that C(t n ) is 
very close to the ground state. 



Modified Langevin Equation 

The ability of this dynamics to maintain the system 
on the BO surface may vary considerably. It is essen- 
tially ideal in systems like water, but potentially some- 
what less perfect in liquid Si at high temperature, where 
swift bonding and rebonding processes take continously 
place. However, in all cases the dynamics is dissipative, 
most likely because the employed propagation scheme 
is not symplectic. Nevertheless, it is possible to rigor- 
ously remedy this downward drift if we assume that the 
forces arising from our dynamics Fpc can be modeled as 
Fpc = Fbo — 7_dR/, which, as we shall see immediately, 
is an excellent assumption. The value of the intrinsic fric- 
tion coefficient does not need to be known but it can 
be bootstrapped by taking a cue from the work of Kra- 
jewski and Parrinello The canonical distribution is 
sampled by using the following Langevin-type equation 



A/ / R / = F pc -7lR-/ + 3/, 



(44) 



where Mj is the ionic mass, jl is a Langevin friction 
coefficient and Ej = af + Ej an additive white noise. 
Using the above assumption Eq. (|4"4f is identically to: 



MjRj = F BO - (jD + 7x) R-i + 3/ 



(45) 



in 



In order to guarantee an accurate sampling of the 
Boltzmann distribution, the noise has to obey the fluc- 
tuation dissipation theorem: 



(0) H 7 (t)) = 6 ( 7J3 + 1L ) M T k B T6 (t) 



(46) 



The choice of jl is arbitrary, while the unknown 
has to be determined by requiring that the aggregate 
noise term generate the correct average temperature, i.e. 
fullfills the equipartition theorem ^|M/R|^ = |fcflT. 

As we see in a moment, this leads to correct a sampling 
of the Boltzmann distribution. In addition, since the 
initial dynamics is quite accurate, 7^ is rather small and 
even dynamical properties can be very well reproduced. 



Illustrative Examples: Liquid Silicon, Silica and 
Water 



For the purpose of demonstrating this new approach, it 
has been implemented in the mixed Gaussian Plane Wave 
(GPW) [zl code Quickstep (z3,1z3|, which ispart of the 
publicly available suite of programs CP2K 75j . In order 
to illustrate that this method works well irrespective of 
band-gap, system size and type, calculations on metallic 
liquid silicon and liquid silica are presented. Both sys- 
tems are known to be very difficult, and are examples of 
liquid metals (Si) as well as of complex, highly polariz- 
able, ionic liquids (SiO^). Furthermore, the simulations 
have been performed at 3000 K and 3500 K respectively, 
which leads to rapidly varying density matrix elements, 
thus making the propagation of the electronic degrees 
of freedom particularly challenging. Hence, the selected 
test cases can be considered as worst-case scenarios for 
any method. 

All simulations have been performed at their exper- 
imental liquid densities using double-zeta valence po- 
larization (DZVP) basis sets, adequate density cut- 



offs, Goedecker-Teter-Hutter pseudopotentials 76-78J 
and the local-density approximation to the exact ex- 
change and correlation functional. For simplicity the 
Brillouin zone is sampled at the T-point only, while 
Eq. (|4"S"]) is integrated using the algorithm of Ricci and 
Ciccotti [z||, with a time step of ft. = 1.0 fs. The friction 
coefficient jl was set equal to zero, while the values for 
7d turned out to be in the range of 10~ 4 fs -1 . The new 
C's are predicted using K = 4 in Eq. (j39|) . which ensures 
time-reversibility up to 0(h e ). 

First, the accuracy in terms of the energetical devia- 
tion from the BO surface is considered. As can be seen 
in FIG. [T] the energies are an upper bound to the ground 
state and are displaced by a very small and approxi- 
mately constant amount. It is also shown that, as al- 
ready mentioned, the deviation from the BO surface can 
be even further reduced by increasing the number of cor- 
rector steps. In fact, it is actually possible to control the 
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FIG. 1. Deviations from the BO surface of liquid SiC>2 with 
respect to total energies (upper panel) and mean force devia- 
tions (lower panel). The deviation in the energies corresponds 
to a constant shift of 4.16 x 10 -4 Hartree per atom for one cor- 
rector step and 3.5 x 10 -5 Hartree per atom for two corrector 
steps. The average mean force deviation is unbiased. 
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FIG. 2. Partial pair-correlation functions g(r) of liquid Si 
(upper left panel) and liquid Si02 at 3000 K and 3500 K 
respectively, using a DZVP Gaussian basis set. 



deviation from the BO surface by varying the number of 
corrector steps in order to achieve a preassigned accuracy 
level. However, in the following only simulations based 
on a single corrector step will be reported. 

Nevertheless, let us now to turn to more realistic prob- 
lems such as those shown in FIG. [2] Although these sim- 
ulations have been performed with only a single corrector 
step, they are still amazingly close to the BOMD refer- 
ence results. It should be emphasized that even in liquid 
Si, which is metallic and poses problems when using an 
ordinary CP scheme, a single corrector step is sufficient. 
This establishes the efficiency of this method, since only 
a single preconditioned gradient calculation with no ad- 
ditional minimization step has to be performed. The 
possible acceleration, in comparison with regular BOMD 
calculations, depends crucially on the system studied. In 
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FIG. 3. The kinetic energy distribution calculated from a 
1 ns trajectory with a time step of 3.25 fs of metallic liquid 
Si64 using a DZVP Gaussian basis set and a density cutoff of 
100 Ry (left panel). Velocity autocorrelation function (upper 
right) and its Fourier transform (lower right) of 32 Water at 
325 K using a TZV2P Gaussian basis set, a density cutoff of 
280 Ry and the BLYP exchange-correlation functional. The 
Langevin friction coefficients are — and jd ~ 10~ 8 fs -1 . 



the undoubtedly difficult cases just presented a speed- 
up of two orders of magnitude compared to using a pure 
extrapolation scheme have been observed. For simpler 
problems still an increase in efficiency of at least one or- 
der of magnitude can be expected. 

In FIG. [3] displays results, which prove that also dy- 
namical properties can be evaluated with accuracy. To 
that extend the velocity autocorrelation function and its 
Fourier transform at 325 K is presented. The results are 
in good agreement with accurate reference calculations 
and are consistent with experiment, as well as ab-initio 
all-electron calculations [SOj , showing that in spite of the 
stochastic nature of Eq. (j45|) dynamical properties can 
also be simulated. This implies, that also chemical reac- 
tions and even non-equilibrium processes can be treated. 
In the same picture it is explicitly verified that the pre- 
vious assumptions are justified, and indeed a canonical 
sampling is performed, by showing that the kinetic en- 
ergy distribution is Maxwellian distributed. To this end, 
a 64 atom liquid Si simulation is carried out for as long as 
1 ns, to reduce the noise and to ensure a proper sampling 
of the relevant kinetic energy distribution tails. 

Due to space considerations only a fraction of the sys- 
tems studied is reported here. Nevertheless, in all cases 
this method has proven to be accurate and the gain in 
speed has always been remarkable 81-86{. Structure re- 
laxation via dynamic annealing and geometry optimiza- 
tion have also been successfully performed [87H92J . Con- 
trary to CPMD and related methods integration time 
steps up to the ionic resonance can be used. Thanks to 
this development it is now possible to perform AIMD sim- 
ulations on medium-sized systems up to a few nanosec- 
onds, thus making a new class of problems accessible to 



ab-initio simulations. 



CONCLUSION 

To conclude it should be noted, that with increasing 
length and time scales CPMD-based approaches are ex- 
pected to become more advantageous than BOMD, since 
otherwise meeting the more and more stringent accu- 
racy requirements of longer simulations and larger system 
sizes would entail an ever tighter wavefunction conver- 
gence. The Car-Parrinello-like approach to BOMD [(3(| 
just described extends the scope of ab-initio simulations 
by combining the best of either method and allows for 
AIMD simulations previously thought not feasible. 
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